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Abstract 

We present a spectral weight conserving formalism for Fermionic thermal Green's functions that 
are discretized in imaginary time r and thus periodic in imaginary ( "Matsubara" ) frequency iuUn- 
The formalism requires a generalization of the Dyson equation G(Go, ^) and the Baym-Kadanoff- 
Luttinger-Ward functional for the free energy (3 ft — r(G). A conformal transformation is used 
to analytically continue the periodized Matsubara Green's function to real frequencies in a way 
that conserves the discontinuity at t = of the corresponding real-time Green's function. This 
allows numerical Green's function calculations of very high precision and it appears to give a well 
controlled convergent approximation in the r discretization. The formalism is tested on dynamical 
mean field theory calculations of the paramagnetic Hubbard model. 



PACS numbers: 



The analytic properties of finite temperature Green's functions is a cornerstone of the 
theory of quantum many body theory. [H [2] The Green's functions can be represented either 
in continuous complex time or by discrete values on an infinite set of Matsubara frequencies. 
The physical retarded Green's function and the corresponding spectral function is given by 
the analytic continuation of this discrete function of Matsubara frequencies to a continuous 
function on the real frequency axis. Although it is extremely elegant the method poses a 
numerical challenge. It was early recognized that the Pade series could be used fit data from a 
finite number of Matsubara frequencies. [3j This scheme is numerically quite ill conditioned 
and much eflFort has been devoted to resolving this difficulty. For quantum Monte Carlo 
methods the Maximum Entropy method [4j that directly computes the spectral function as 
a probabilility distribution is instead widely used. 

We develop a method to work with a Fermion Green's function defined exactly only at 

equally spaced points tj = j^j in imaginary time. In this space of dimension N a discrete 
solution is computed numerically exactly. A conformal transformation is used to construct 
an exact analytic continuation using a rational function. Since N does not have to be large 
to yield stable results, calculations can be done using very high precision. This resolves 
many of the difficulties others have had with an ill conditioned Pade series. [5j 

Our approach assures that several basic conditions are satisfied. In addition to the obvious 
demand that the limit A ^ oc yields a proper continuum limit, we obtain three additional 
properties valid for all A that are not present in previous approaches and which give a 
well controlled and rapidly convergent result as A is increased. For all values of A we 
find that (a) the Green's function obeys exactly a Luttinger-Ward variational principle 
(b) the free energy is exact for noninteracting particles (c) a numerically exact analytic 
continuation exists that reproduces the data on the imaginary frequency axis and has the 
proper discontinuity and analytic structure at t = 0. 

Let us consider Fermions described by a Hamiltonian H = Hq + V with Hq = ^kC^Ck^ 
and where c|. represents Fermion creation operators of a state \k) and where k may 
be momentum and spin. The chemical potential /x is absorbed in Ck- The interac- 
tion 18 V = ^^^kk' qq'^kk'qq'clc\,CqfCq. Wc cousldcr the ouc-particle Green's function 
Gk{r,r') = -{T{ck{r)cl{r')) = -^Tr{e-^^T{e^^ CkC'^^ e^'^ cle'^'^)} where Tr is the sum 
over a complete set of states, T is time ordering, (5 is inverse temperature, and Z is the par- 
tition function Z = Tr{e~^^}. We assume G is diagonal in k and there are no anomalous 
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''superconducting" terms. The function Gk{T^T') is defined on —(5 < r — r' < (3 and obeys 
Gk{(^ — r) = —Gk{—T). This antiperiodicity condition allows Gkir) to be transformed to 
Matsubara frequencies iujn = + |) for integer n. 

The expansion of Gkir) in powers of V can be expressed as sums of connected diagrams 
consisting of the vertex at a time r and the non-interacting Green's function Go^k{^ — ^^) = 
— {T{ck{T)c\T'))Q^ where internal times are integrated over as dr and an n'th order 
diagram has a prefactor ^ with ni the number of Fermion loops. [2j In this work we 

discretize r and replace integrals over continuous time < r < ;5 with sums over r^, by 
carefully defining a discretized Green's function and self energy. 

We start by defining V = ^ = ^ and discretize the non-interacting Green's 

function Go,fc(r) = e— ^[n;(e,)e(-r+0+) + (n;(e,)-l)e(r-0+)], with n/(e,) = l/(e^^^ + l), 
to times tj = 2r]j . The value at tq = is defined as the average of the two limits so 
that Go,/c(ro) = nf{ek) - \. We can expand G^^^{t^) = ^Y.n=o ^''"^"^'^oA^^n) and the 
''Matsubara transform" G{i(jJn) = ^{G{{rj})) is 

N-l 

Go^iuJn) = ^vYl ^'""^^'Go^kiTj) = 77coth r]{iujn - Ck) . (1) 
j=o 

For compactness we will usually drop the index k and frequency iuJn when it is clear from 
the context. We see that Go is periodic [6j under iuj^ i(^n + i^N and in the limit ^ oc 
{rj 0) the continuum expression Go = — ^k) is appropriately recovered. We also 

define the additional Green's functions G^ = GQ±r] in analogy to the r = 0^ limit of Gq{t). 

Let us now define the periodized full Green's function Gk{ic^n) and the self-energy T>k{i^n) 
through the two expressions 

Gk{iuJn) = r]cothr]{iu)n - - T>k{i^n)) (2) 

and T>k{i^n) = sGt{iuj ) ' object $ is the functionaljll l8j| defined as the sum of linked 
closed skeleton diagrams of Gk{ioOn)-> except for the 1st order diagrams where we use 
G^{i(jUn) = Gk{i(jOn) + T]. The self energy is thus given by the amputated skeleton diagrams 
and the resulting set of equations using Eq. [2j 

Equation [2] is a generalization of the standard Dyson equation G^^{i(jJn) = iuOn — — 
^k{i^n)) and reduces to the latter as ^ 0. It also preserves the property that a constant, k 
and ujn independent, self energy acts a chemical potential. Defining S = - tanh {r]T>k{i(jOn)) 
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which reduces to Tik{iuJn) in the hmit ry ^ 0. Eq. [2| can also be rewritten in the more 
suggestive form 



= . , (3) 



Go' - S 
1- 772^0 

which again reduces to the ordinary Dyson equation for 77 ^ 0. 

Consider now the free energy Q = —^InZ. As shown by Luttinger and Ward [7], for 
the standard formahsm, the free energy can be expressed as /3Q = F = $({G}) — TrT>G + 
Tr log(—G), where Tr = J^kun' expression F = F({G}) also provides a variational for- 
mulation where F = /3Q corresponds to a stationary point 6r/6G = that yields the Dyson 
equation = Gq^ — S. We now describe how to generalize this variational formulation 
to be consistent with Eq. [sjand demands (a)-(b) stated in the introduction. 

We define the following generalization 

F = $({G}) - rr(G+S) + Tr log {-G-/{2r])) (4) 

where G^ = G ±r]. The expression F reduces to the standard expression in the limit rj ^ 
and the stationarity condition 5T/5G = gives the generalized Dyson equation Eq. [3j 
In addition, it can be shown that in the noninteracting limit when S = $ = the free 
energy in Eq,4 is exact for all values of through the expression ^^^In ^(— Gq ^(icj^)) = 
— J2k ^^(^ + e~^^^) = /3f]o- The functional F({G}) can also be shown to give a value for the 
total particle number which is consistent with the formalism i.e. ^ — — ^Tr G^ . Baym 
and Kadanofl[ [8j noted that conserving approximations can be found by including only a 
finite number of diagrams in the Luttinger- Ward function. This construction can now be 
used for the discretized Green's functions to include only a subset of diagrams in $. 

Let us now explore the analytic structure of the periodized Green's function. Starting 
with the spectral representation in terms of a complete set of eigenstates H\n) = En\n)^ for 
r > 0, Gk{r) = — ^ Z]rnn^"^^''^^^^''"^'^^K^l4l^)P using Eq. m we can compute G as 



— A{k, uj)r] coth r]{iujn - uj) , (5) 
-00 27r 

where the spectral function is given by the conventional expression A{k^uj) = 
\ Yjmn \{m\cl\n)\'^e~^^''{l + e~f^'^)27r5{Em - - uj). Eq. [5] gives the analytic continu- 
ation io z = UJ + ilm{z) by letting iuj^ z and shows that Gk{z) is analytic except where 
Im{z) is an integer multiple oi Vt^ = ^. Using r]Imcothr]{(jUo + iO^ — u)) = —7r5{uJo — u)) 
and defining G{z = uj + iO^) = G^{uj) and G{z = cj + zO") = G^{uj)^ we formally recover the 



standard expression A{k^(jj) = ±2ImG^^^{(jj). Inserting into Eq. [5] gives the generalized 
Kramers-Kronig relation 

We will now show how to use a conformal transformation to make an analytic continuation 
of the periodized Green's function to a rational function with simple poles. Consider a 
spectral function of the form 

L.Juj) = ^- — — — — — . (7) 

smh r]\uj — 6 + ^7) smhr]\uj — e — i^) 

which reduces to the Lorentzian ^^_'^y^^2 in the limit rj ^ 0. Because of the antiperiodicity 
and to ensure positive spectral weight we can assume < ^ < We can evaluate the 

integral in Eq. [5] by a closed contour containing three poles, as shown in Figure [T] The 
result, for < Imz < Qn^ is 

1 f dz' 

G{z) = -f —L,,,{z')vcothv{z - z') = (8) 

77/77 T] \ 

- f coth -{z — e + i^) + tanh -{z — e — 



and using the same count our integral it can also shown be that L[uj) is properly normalized: 
/-^ ^I/(cl;) = 1. The function G{z) is thus completely free from singularities in the strip 
< Imz < and has poles outside the strip at e — ^7 and e + + i^N and obeys 
G{z) = G{z + 2iQ]\i) as shown in Figure [3^. In the limit 77 ^ the function reduces to 
G^(^) = z-(e-i-f) which is the usual retarded Green's function from a Lorentzian spectral 
weight. Evaluating G{z) in the strip —Qn < I'^z < gives an expression which is analytic 
in that strip and which analogously corresponds to the advanced Green's function. The 
full analytic Green's function with periodically repeated branch cuts corresponds to gluing 
together the two branches as indicated in Figure [3|3. In the limit r] ^ 0^ = tt/t] ^ 00 
thus reproducing the standard structure of the analytic Green's function with a branch cut 
on the real axis. 

Let us now consider a spectral function defined by a set of such "periodized Lorentzians" 
L^^^^^{uj). For < Imz < Qn we find 

G{z) = ^^K coth |(z - 6, + ij,) + (9) 

77 

a* tanh -(z - - ^7^)] , 
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FIG. 1: (Color online) Contour for the integration of Eq. [s] 

where we have allowed for a complex prefactor that must satisfy | ^j^ictu + ctl) = 1 to 
conserve total spectral weight. Given a set of values of Gn = Gk{i(jOn) we would like to 
extract the values of a^^, and 7^ that solves Eq. [9| The latter can be cast into the form 
of a rational function containing only simple poles (as opposed to periodically repeated) by 
means of the conformal mapping 



z coth -z 
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which gives 
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with = coth(^(e^^ — i^j^) — if). The transformation maps the strip < Imz < 2Qn 
the entire complex plane. The strip J^at < Im z < 2Qn that contains the singularities is 
mapped to the interior of the unit circle and the strip < Imz < that is free from 
singularities is outside, as exemplified in Figure ^ and d. The points z = ±oc on the real 
axis map to = ±1 and the point z = iuj(^N -i)/2 maps to = oo. 
The function G{z') should obey a number of properties. Since Eq. 
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can be written as 



a rational function G{z') = const + P(z')/Q{z') where P is an M — 1 degree and Q is an 
M'th degree polynomial of z' and const = — f '^{O'l'Pi' + (i-IpI), we can identify pu with the 
roots of Q and |(1 — pDa^, as the residues oi P/Q. Taking odd gives a precise boundary 
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condition G{z' oo) const G(^n-i)/2- We also require G{z' = 1) = —G{z' = — 1) = 
2 + ^y) = to obtain the proper normalization. A crucial observation is that Eq. 

3] preserves this condition independent of S and demonstrates that total spectral weight is 
preserved by the periodized Dyson equation. By straightforward calculation from Eq. [5] 
we also conclude that G'{z = 1) = G\z = —1) = which further constrains the analytic 
continuation. We thus find P{z^) and Q{z') by fitting G{z') — const to (M — 1) /M Pade form, 
yielding M = (A^ + 3)/2. In addition we have the symmetry G{i(jJN-i-n) = G''{i(jJn) resulting 



in M/2 independent poles. Thus n = (A^ + 3)/4 in Eq. 11, Having identified the parameters 

and in the fit, the spectral function can be evaluated through A{(jj) = —2ImG{uj) using 
Eq|9| For every value of we can also obtain a properly normalized spectral weight using 
A{uj) = 2 ^^^''(^^Ie^2l^j2 ~ ^ resulting in Green's function with the proper discontinuity 
Gk{T = 0-) - Gk{T = 0+) = t^ik, io) = 1. 

We tested our method on the Dynamical mean field theory (DMFT) method using the 
IPT approximation with the half filled paramagnetic Hubbard model. [914T2] We found that 
the convergence was significantly enhanced by evaluating the dynamic part of the self energy 
as tdLiih r]T>{iuJn) = vU^^ {Gl^impi^j)) where Eq. [ijhas been used. This expression ensures 
that /mS lies within the appropriate bounds and by expanding tanh in powers of S it 
is equivalent to the standard expression to order [/^. The detailed implications of this 
procedure should be explored further. 

Figure [2] shows the spectral function for varying bare energy A{ek^uj) at ;5 = 25 (T = 
VI^/50), in the metallic, U = 2 (with N = 25), and insulating, U = A {N = 45) phases using 
the standard semicircular bare density of states[12j and 40 significant digits to compute the 
rational function. This gives 7 and 12 independent complex poles respectively. The location 
of poles are indicated in Fig. [2] for Ck = 0. Since S is fc-independent we can extract it with a 
single Fade fit and use this to generate the spectral function for any Ck- We have found that 
the number of frequencies needed is roughly A > /5 for convergence of the IPT recursion 
algorithm and as further test we have used the method at low temperature /3 > 500 and a 
number of significant digits roughly also equal to A. 

As a second test of the analytic continuation method, we considered a noninteracting 
single impurity Anderson model whose spectral function consists of a sharp resonance as 
well as a continuum. These features have been found to be difficult to reproduce in detail 
using ordinary Fade methods. [13] We used Eq. [5] to compute Gimp{i(^n) from the exact 
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FIG. 2: (Color online) Spectral functions A{ek,(jj) at /3 = 25 for metal ?7 = 2 (a) and insulator 
(7 = 4 (b) for bare energies Ck = —1, 0, 1. (Dashed curves are /3 = 500 and /3 = 200 respectively.) 
The corresponding Im{Ge^=o{z')) is shown in (c) and (d) including circles marking the location of 
poles. The poles are located inside the unit circle in accordance with Eq. [TT] Empty dots represent 
a subset of the Matsubara values {z'{iujn)). The color scale is light for large positive and negative 
values with dark colors near zero. 

^impi"^) using = 61, /3 = 25, and 60 significant figures. Using this we find 16 independent 
complex poles that to the eye reproduces the exact spectral function. It has an integrated 
rms deviation over total weight (excluding the resonance) of 4 • 10~^. The resonance has 
a width of 10~^ and a normalization which is within 10~^ of the exact value. A major 
reason for our success is the high precision for the Fade fit rather than a large number 
of poles. Fitting a greater number of less accurate data points does not yield comparable 
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accuracy. This motivates using our method to compute a Green's function to extremely high 
precision in a relatively small dimensional space and to use the combination of conformal 
transformation and Pade fit to infer the analytic continuation to the real axis. 

We have presented a method to work with discrete-time Matsubara Green's functions. 
In spite of working in a finite dimensional space the method yields an analytic continuation 
that obeys the boundary condition at r = 0, a problem that has plagued other previous 
discretization schemes. Within this small space, numerical calculations can be done to the 
extremely high accuracy that is necessary to numerically obtain a meaningful analytic con- 
tinuations from imaginary to the real time axis. The method should have wide applicability 
to all problems requiring numerical evaluation of Matsubara Green's function in condensed 
matter physics nuclear physics and quantum field theory. 
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